{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.io\n",
    "from scipy.stats import itemfreq\n",
    "from statistics import mode\n",
    "import time\n",
    "import multiprocessing as mp\n",
    "\n",
    "import pydot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decision Tree class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class decision_tree:\n",
    "    \n",
    "    class node:\n",
    "        def __init__(self, left, right, split_rule, is_leaf, label):\n",
    "            self.left = left\n",
    "            self.right = right\n",
    "            self.split_rule = split_rule\n",
    "            self.is_leaf = is_leaf\n",
    "            self.label = label\n",
    "    \n",
    "    def __init__(self, max_depth=1e10):\n",
    "        self.max_depth = max_depth\n",
    "\n",
    "    def max_count(self, array):\n",
    "        return max(set(array), key=array.tolist().count)\n",
    "        \n",
    "    # utility function for entropy calculation\n",
    "    def entropy(self, indices):\n",
    "        p = itemfreq(a)[:, 1].astype(float) / len(indices)\n",
    "        return -p.dot(np.log2(p))\n",
    "    \n",
    "    # calculate entropy the number of instances in each class in known\n",
    "    def entropy_n(self, all_n):\n",
    "        p = all_n / (np.sum(all_n)+1e-20)\n",
    "        return -p.dot(np.log2(p+1e-20))\n",
    "    \n",
    "    # calculate the impurity(\"badness\") of the specified split on the input data\n",
    "    def impurity(self, left_label_hist, right_label_hist):\n",
    "        Sl = np.sum(left_label_hist)\n",
    "        Sr = np.sum(right_label_hist)\n",
    "        return (Sl*self.entropy_n(left_label_hist) + Sr * self.entropy_n(right_label_hist)) / (Sl+Sr)\n",
    "    \n",
    "    # find the threshold that best split the data points with a certain feature\n",
    "    # Note: <= th goes to S_left and > th goes to S_right\n",
    "    def find_threshold(self, feature, labels):\n",
    "        all_f = sorted(set(feature)) # sorted in ascending order\n",
    "        all_l = set(labels) # list unique labels\n",
    "        \n",
    "        freq_mat = np.zeros([len(all_f), len(all_l)])\n",
    "        for i, f in enumerate(all_f):\n",
    "            for j, l in enumerate(all_l):\n",
    "                freq_mat[i, j] = len(labels[np.where(labels[np.where(feature==f)]==l)])\n",
    "        \n",
    "        # calculate the average of two neighboring values as threshold\n",
    "        # iterates from min to max\n",
    "        all_threshold = (np.hstack((all_f[1:], all_f[-1])) + all_f) / 2.\n",
    "        \n",
    "        # in the beginning, all goes to the right node\n",
    "        n_left = np.zeros([len(all_l)])\n",
    "        n_right = np.sum(freq_mat, axis=0)\n",
    "        n_left_sum = 0\n",
    "        min_threshold = all_threshold[0]\n",
    "        min_H = self.impurity(n_left, n_right)\n",
    "        # loop through all threshold to find the one with the minimum impurity\n",
    "        for i, th in enumerate(all_threshold):\n",
    "            n_left += freq_mat[i, :]\n",
    "            n_right -= freq_mat[i, :]\n",
    "            H = self.impurity(n_left, n_right)\n",
    "            if H < min_H:\n",
    "                min_H = H\n",
    "                min_threshold = th\n",
    "        return min_threshold, min_H\n",
    "    \n",
    "    \n",
    "    # find the best feature and threshold to split data points\n",
    "    def segmenter(self, data, labels, m=-1):\n",
    "        d = data.shape[1]\n",
    "        if m == -1:\n",
    "            all_features = np.arange(d)\n",
    "        else:\n",
    "            all_features = np.random.choice(range(d), m, replace=False)\n",
    "        min_H = 1e20\n",
    "        min_th = 0\n",
    "        min_i = 0\n",
    "        for i in all_features:\n",
    "            threshold, H = self.find_threshold(data[:, i], labels)\n",
    "            if H < min_H:\n",
    "                min_H = H\n",
    "                min_th = threshold\n",
    "                min_i = i\n",
    "        return min_i, min_th\n",
    "    \n",
    "    # the recurrence function that builds the decision tree\n",
    "    def grow_tree(self, S, depth, m=-1):\n",
    "        if len(set(self.labels[S])) == 1 or depth >= self.max_depth: # pure node or reach maximum depth\n",
    "            return self.node(left=None, right=None, split_rule=None, is_leaf=1, \\\n",
    "                             label=self.max_count(self.labels[S]))\n",
    "        else:\n",
    "            min_i, min_th = self.segmenter(self.data[S, :], self.labels[S], m=m)\n",
    "            Sl = [j for j in S if self.data[j, min_i] <= min_th]\n",
    "            Sr = [j for j in S if self.data[j, min_i] > min_th]\n",
    "            if len(Sl) == 0 or len(Sr) == 0:\n",
    "                return self.node(left=None, right=None, split_rule=None, is_leaf=1, \\\n",
    "                                 label=self.max_count(self.labels[S]))\n",
    "            else:\n",
    "                return self.node(left=self.grow_tree(Sl, depth+1, m), right=self.grow_tree(Sr,depth+1, m), \\\n",
    "                                 split_rule = (min_i, min_th), \\\n",
    "                            is_leaf=0, label=None)\n",
    "        \n",
    "    # train the decision tree\n",
    "    def train(self, data, labels, m=-1):\n",
    "        self.data = data\n",
    "        self.labels = labels\n",
    "        S = np.array(range(len(labels)))\n",
    "        self.root = self.grow_tree(S, 1, m=m)\n",
    "        self.data = None\n",
    "        self.labels = None\n",
    "        return self\n",
    "        \n",
    "    # predict labels of test data\n",
    "    def predict(self, data, verbose=0):\n",
    "        if data.ndim == 1: # special case of only 1 row (it becomes a 1d vector in numpy)\n",
    "            data = np.reshape(data, [1, len(data)])\n",
    "            N = 1\n",
    "        else:\n",
    "            N = data.shape[0]\n",
    "        labels = np.zeros(N)    \n",
    "            \n",
    "        # predict each data point    \n",
    "        for i in range(N):\n",
    "            d = data[i, :]\n",
    "            current_node = self.root\n",
    "            depth = 0\n",
    "            # going down along the tree\n",
    "            while not current_node.is_leaf: # not reach leaf yet\n",
    "                idx = current_node.split_rule[0]\n",
    "                th = current_node.split_rule[1]\n",
    "                if d[idx] <= th:\n",
    "                    current_node = current_node.left\n",
    "                    if verbose == 1 and depth == 0:\n",
    "                        print('Feature %d <= %f' % (idx, th))\n",
    "                    if verbose == 2:\n",
    "                        print('Feature %d <= %f' % (idx, th))\n",
    "                else:\n",
    "                    current_node = current_node.right\n",
    "                    if verbose == 1 and depth == 0:\n",
    "                        print('Feature %d > %f' % (idx, th))\n",
    "                    if verbose == 2:\n",
    "                        print('Feature %d > %f' % (idx, th))\n",
    "                depth += 1\n",
    "            if verbose == 2:\n",
    "                print('Classification: %d' % current_node.label)\n",
    "                \n",
    "            labels[i] = current_node.label\n",
    "        return labels\n",
    "\n",
    "    # calculate the prediction accuracy\n",
    "    def accuracy(self, data, true_labels):\n",
    "        labels = self.predict(data, verbose=False)\n",
    "        N = len(labels)\n",
    "        return np.sum(labels == true_labels) / float(N)\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Forest Class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# use multiproccessing to speed up\n",
    "class random_forest:\n",
    "    def __init__(self, n_trees=20, n_sample=1000, n_feature=-1, max_depth=1e10):\n",
    "        self.n_trees = n_trees\n",
    "        self.n_sample = n_sample\n",
    "        self.n_feature = -1\n",
    "        self.max_depth = max_depth\n",
    "        self.trees = np.array([decision_tree(max_depth)] * n_trees)\n",
    "        \n",
    "    def train_parallel(self, data, labels, verbose=True):\n",
    "        if self.n_feature == -1:\n",
    "            d = data.shape[1]\n",
    "            self.n_feature = int(np.sqrt(d)) # num of random features = sqrt(d) is a good guess\n",
    "        pool = mp.Pool()\n",
    "        results = np.zeros(self.n_trees, dtype=object)\n",
    "        print('# Trees = %d' % self.n_trees)\n",
    "        for i, dt in enumerate(self.trees):\n",
    "#             print(i)\n",
    "            idx = np.random.choice(range(len(data)), self.n_sample)\n",
    "            sub_data = data[idx, :]\n",
    "            sub_labels = labels[idx]\n",
    "            results[i] = pool.apply_async(dt.train, args=(sub_data, sub_labels, self.n_feature,))\n",
    "            if verbose:\n",
    "                print('%d: Tree loaded.' % i)\n",
    "#             dt.train(sub_data, sub_labels, m=self.n_feature) # activate the random feature mode\n",
    "        self.trees = [p.get() for p in results]\n",
    "        print('Trees trained.')\n",
    "        pool.close()\n",
    "        \n",
    "    def train(self, data, labels):\n",
    "        if self.n_feature == -1:\n",
    "            d = data.shape[1]\n",
    "            self.n_feature = int(np.sqrt(d)) # num of random features = sqrt(d) is a good guess\n",
    "        print('# Trees = %d' % self.n_trees)\n",
    "        for i, dt in enumerate(self.trees):\n",
    "            print(i)\n",
    "            idx = np.random.choice(range(len(data)), self.n_sample)\n",
    "            sub_data = data[idx, :]\n",
    "            sub_labels = labels[idx]\n",
    "            dt.train(sub_data, sub_labels, m=self.n_feature) # activate the random feature mode\n",
    "    \n",
    "    def predict(self, data, verbose=0):\n",
    "        if data.ndim == 1: # special case of only 1 row (it becomes a 1d vector in numpy)\n",
    "            data = np.reshape(data, [1, len(data)])\n",
    "            N = 1\n",
    "        else:\n",
    "            N = data.shape[0]\n",
    "        labels = np.zeros(N)    \n",
    "            \n",
    "        # predict each data point    \n",
    "        for i in range(N):\n",
    "            votes = np.zeros(self.n_trees)\n",
    "            for j, t in enumerate(self.trees):\n",
    "                votes[j] = t.predict(data[i, :], verbose=verbose)\n",
    "            labels[i] = t.max_count(votes)\n",
    "        return labels\n",
    "    \n",
    "    def accuracy(self, data, true_labels):\n",
    "        labels = self.predict(data, verbose=False)\n",
    "        N = len(labels)\n",
    "        return np.sum(labels == true_labels) / float(N)\n",
    "     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Utility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def normalize_row(X):\n",
    "    Xn = np.zeros(X.shape)\n",
    "    for i in range(X.shape[0]):\n",
    "        x = X[i, :]\n",
    "        Xn[i, :] = (x/(np.sqrt(x.dot(x))+1e-15))\n",
    "    return Xn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def max_count(lst):\n",
    "    return max(set(lst), key=lst.tolist().count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. spam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "23702"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = scipy.io.loadmat('spam_dist/spam_data.mat')\n",
    "train_X = data['training_data']\n",
    "train_y = data['training_labels'].ravel()\n",
    "test_X = data['test_data']\n",
    "len(train_X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(18961, 32)\n"
     ]
    }
   ],
   "source": [
    "S = np.random.choice(range(len(train_y)), int(len(train_y)/5*4), replace=False)\n",
    "train_data = train_X[S, :]\n",
    "train_labels = train_y[S]\n",
    "S = np.setdiff1d(range(len(train_y)), S)\n",
    "validation_data = train_X[S, :]\n",
    "validation_labels = train_y[S]\n",
    "print(train_data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Decision Tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total training time:  9.92026662826538\n",
      "Accuracy:  0.806866726438\n",
      "Accuracy:  0.793714406243\n"
     ]
    }
   ],
   "source": [
    "dt = decision_tree(15)\n",
    "start = time.time()\n",
    "dt.train(train_data, train_labels)\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)\n",
    "\n",
    "\n",
    "l = dt.accuracy(train_data, train_labels)\n",
    "print('Accuracy: ', l)\n",
    "\n",
    "l = dt.accuracy(validation_data, validation_labels)\n",
    "print('Accuracy: ', l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feature 28 <= 0.500000\n",
      "Feature 29 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 25 <= 1.500000\n",
      "Feature 7 <= 0.500000\n",
      "Feature 13 <= 0.500000\n",
      "Feature 3 <= 0.500000\n",
      "Feature 0 <= 0.500000\n",
      "Feature 26 <= 0.500000\n",
      "Feature 31 <= 0.500000\n",
      "Feature 6 <= 0.500000\n",
      "Feature 20 <= 0.500000\n",
      "Feature 16 <= 0.500000\n",
      "Feature 30 <= 0.500000\n",
      "Classification: 0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.])"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dt.predict(train_X[0, :], verbose=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Random Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Trees = 20\n",
      "0: Tree loaded.\n",
      "1: Tree loaded.\n",
      "2: Tree loaded.\n",
      "3: Tree loaded.\n",
      "4: Tree loaded.\n",
      "5: Tree loaded.\n",
      "6: Tree loaded.\n",
      "7: Tree loaded.\n",
      "8: Tree loaded.\n",
      "9: Tree loaded.\n",
      "10: Tree loaded.\n",
      "11: Tree loaded.\n",
      "12: Tree loaded.\n",
      "13: Tree loaded.\n",
      "14: Tree loaded.\n",
      "15: Tree loaded.\n",
      "16: Tree loaded.\n",
      "17: Tree loaded.\n",
      "18: Tree loaded.\n",
      "19: Tree loaded.\n",
      "Trees trained.\n",
      "Total training time:  5.115397930145264\n",
      "Accuracy:  0.814250303254\n",
      "Accuracy:  0.810166631512\n"
     ]
    }
   ],
   "source": [
    "n_trees=20\n",
    "n_sample=18000\n",
    "n_feature=-1\n",
    "max_depth=30\n",
    "rf = random_forest(n_trees=n_trees, n_sample=n_sample, n_feature=n_feature, max_depth=max_depth)\n",
    "start = time.time()\n",
    "rf.train_parallel(train_data, train_labels)\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)\n",
    "\n",
    "l = rf.accuracy(train_data, train_labels)\n",
    "print('Accuracy: ', l)\n",
    "\n",
    "l = rf.accuracy(validation_data, validation_labels)\n",
    "print('Accuracy: ', l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 16 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n",
      "Feature 19 <= 0.500000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.])"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf.predict(train_X[0, :], verbose=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Kaggle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer\n",
    "# test tf-idf feature\n",
    "from sklearn.feature_extraction.text import TfidfVectorizer\n",
    "import glob\n",
    "SPAM_DIR = 'spam/'\n",
    "HAM_DIR = 'ham/'\n",
    "TEST_DIR = 'test/'\n",
    "NUM_TEST_EXAMPLES = 10000\n",
    "spam_filenames = glob.glob('spam_dist/' + SPAM_DIR + '*.txt')\n",
    "ham_filenames = glob.glob('spam_dist/' + HAM_DIR + '*.txt')\n",
    "test_filenames = ['spam_dist/' + TEST_DIR + str(x) + '.txt' for x in range(NUM_TEST_EXAMPLES)]\n",
    "\n",
    "all_text = []\n",
    "for file in spam_filenames+ham_filenames: # use only training set data to build BOG\n",
    "    with open(file, \"r\", encoding='utf-8', errors='ignore') as f:\n",
    "        all_text.append(f.read())\n",
    "        \n",
    "all_test_text = []\n",
    "for file in test_filenames: # use only training set data to build BOG\n",
    "    with open(file, \"r\", encoding='utf-8', errors='ignore') as f:\n",
    "        all_test_text.append(f.read())\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# vectorizer = CountVectorizer(min_df=4) # min word length=4\n",
    "vectorizer = TfidfVectorizer()\n",
    "train_X = normalize_row(vectorizer.fit_transform(all_text).toarray())\n",
    "test_X = normalize_row(vectorizer.transform(all_test_text).toarray())\n",
    "train_y = np.concatenate((np.ones(len(spam_filenames)), np.zeros(len(ham_filenames))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(23702, 2000)\n"
     ]
    }
   ],
   "source": [
    "std = np.std(train_X, axis=0)\n",
    "# use the max N variance word to build model\n",
    "n_max = 2000\n",
    "max_idx = std.argsort()[-n_max:]\n",
    "train_X_max = train_X[:, max_idx]\n",
    "test_X_max = test_X[:, max_idx]\n",
    "print(train_X_max.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "S = np.random.choice(range(len(train_y)), int(len(train_y)/5*4), replace=False)\n",
    "train_data = train_X_max[S, :]\n",
    "train_labels = train_y[S]\n",
    "S = np.setdiff1d(range(len(train_y)), S)\n",
    "validation_data = train_X_max[S, :]\n",
    "validation_labels = train_y[S]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Trees = 128\n",
      "0: Tree loaded.\n",
      "1: Tree loaded.\n",
      "2: Tree loaded.\n",
      "3: Tree loaded.\n",
      "4: Tree loaded.\n",
      "5: Tree loaded.\n",
      "6: Tree loaded.\n",
      "7: Tree loaded.\n",
      "8: Tree loaded.\n",
      "9: Tree loaded.\n",
      "10: Tree loaded.\n",
      "11: Tree loaded.\n",
      "12: Tree loaded.\n",
      "13: Tree loaded.\n",
      "14: Tree loaded.\n",
      "15: Tree loaded.\n",
      "16: Tree loaded.\n",
      "17: Tree loaded.\n",
      "18: Tree loaded.\n",
      "19: Tree loaded.\n",
      "20: Tree loaded.\n",
      "21: Tree loaded.\n",
      "22: Tree loaded.\n",
      "23: Tree loaded.\n",
      "24: Tree loaded.\n",
      "25: Tree loaded.\n",
      "26: Tree loaded.\n",
      "27: Tree loaded.\n",
      "28: Tree loaded.\n",
      "29: Tree loaded.\n",
      "30: Tree loaded.\n",
      "31: Tree loaded.\n",
      "32: Tree loaded.\n",
      "33: Tree loaded.\n",
      "34: Tree loaded.\n",
      "35: Tree loaded.\n",
      "36: Tree loaded.\n",
      "37: Tree loaded.\n",
      "38: Tree loaded.\n",
      "39: Tree loaded.\n",
      "40: Tree loaded.\n",
      "41: Tree loaded.\n",
      "42: Tree loaded.\n",
      "43: Tree loaded.\n",
      "44: Tree loaded.\n",
      "45: Tree loaded.\n",
      "46: Tree loaded.\n",
      "47: Tree loaded.\n",
      "48: Tree loaded.\n",
      "49: Tree loaded.\n",
      "50: Tree loaded.\n",
      "51: Tree loaded.\n",
      "52: Tree loaded.\n",
      "53: Tree loaded.\n",
      "54: Tree loaded.\n",
      "55: Tree loaded.\n",
      "56: Tree loaded.\n",
      "57: Tree loaded.\n",
      "58: Tree loaded.\n",
      "59: Tree loaded.\n",
      "60: Tree loaded.\n",
      "61: Tree loaded.\n",
      "62: Tree loaded.\n",
      "63: Tree loaded.\n",
      "64: Tree loaded.\n",
      "65: Tree loaded.\n",
      "66: Tree loaded.\n",
      "67: Tree loaded.\n",
      "68: Tree loaded.\n",
      "69: Tree loaded.\n",
      "70: Tree loaded.\n",
      "71: Tree loaded.\n",
      "72: Tree loaded.\n",
      "73: Tree loaded.\n",
      "74: Tree loaded.\n",
      "75: Tree loaded.\n",
      "76: Tree loaded.\n",
      "77: Tree loaded.\n",
      "78: Tree loaded.\n",
      "79: Tree loaded.\n",
      "80: Tree loaded.\n",
      "81: Tree loaded.\n",
      "82: Tree loaded.\n",
      "83: Tree loaded.\n",
      "84: Tree loaded.\n",
      "85: Tree loaded.\n",
      "86: Tree loaded.\n",
      "87: Tree loaded.\n",
      "88: Tree loaded.\n",
      "89: Tree loaded.\n",
      "90: Tree loaded.\n",
      "91: Tree loaded.\n",
      "92: Tree loaded.\n",
      "93: Tree loaded.\n",
      "94: Tree loaded.\n",
      "95: Tree loaded.\n",
      "96: Tree loaded.\n",
      "97: Tree loaded.\n",
      "98: Tree loaded.\n",
      "99: Tree loaded.\n",
      "100: Tree loaded.\n",
      "101: Tree loaded.\n",
      "102: Tree loaded.\n",
      "103: Tree loaded.\n",
      "104: Tree loaded.\n",
      "105: Tree loaded.\n",
      "106: Tree loaded.\n",
      "107: Tree loaded.\n",
      "108: Tree loaded.\n",
      "109: Tree loaded.\n",
      "110: Tree loaded.\n",
      "111: Tree loaded.\n",
      "112: Tree loaded.\n",
      "113: Tree loaded.\n",
      "114: Tree loaded.\n",
      "115: Tree loaded.\n",
      "116: Tree loaded.\n",
      "117: Tree loaded.\n",
      "118: Tree loaded.\n",
      "119: Tree loaded.\n",
      "120: Tree loaded.\n",
      "121: Tree loaded.\n",
      "122: Tree loaded.\n",
      "123: Tree loaded.\n",
      "124: Tree loaded.\n",
      "125: Tree loaded.\n",
      "126: Tree loaded.\n",
      "127: Tree loaded.\n",
      "Trees trained.\n",
      "Total training time:  213.49836039543152\n"
     ]
    }
   ],
   "source": [
    "# tune random forest parameters\n",
    "n_trees=128 # 128 is good\n",
    "n_sample=18000\n",
    "n_feature=-1 # -1 ==> sqrt(d)\n",
    "max_depth=80 # deeper tree has lower bias\n",
    "\n",
    "rf = random_forest(n_trees=n_trees, n_sample=n_sample, n_feature=n_feature, max_depth=max_depth)\n",
    "\n",
    "start = time.time()\n",
    "rf.train(train_data, train_labels) # use all data\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy:  0.974267032272\n"
     ]
    }
   ],
   "source": [
    "l = rf.accuracy(validation_data, validation_labels)\n",
    "print('Accuracy: ', l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Trees = 128\n",
      "0: Tree loaded.\n",
      "1: Tree loaded.\n",
      "2: Tree loaded.\n",
      "3: Tree loaded.\n",
      "4: Tree loaded.\n",
      "5: Tree loaded.\n",
      "6: Tree loaded.\n",
      "7: Tree loaded.\n",
      "8: Tree loaded.\n",
      "9: Tree loaded.\n",
      "10: Tree loaded.\n",
      "11: Tree loaded.\n",
      "12: Tree loaded.\n",
      "13: Tree loaded.\n",
      "14: Tree loaded.\n",
      "15: Tree loaded.\n",
      "16: Tree loaded.\n",
      "17: Tree loaded.\n",
      "18: Tree loaded.\n",
      "19: Tree loaded.\n",
      "20: Tree loaded.\n",
      "21: Tree loaded.\n",
      "22: Tree loaded.\n",
      "23: Tree loaded.\n",
      "24: Tree loaded.\n",
      "25: Tree loaded.\n",
      "26: Tree loaded.\n",
      "27: Tree loaded.\n",
      "28: Tree loaded.\n",
      "29: Tree loaded.\n",
      "30: Tree loaded.\n",
      "31: Tree loaded.\n",
      "32: Tree loaded.\n",
      "33: Tree loaded.\n",
      "34: Tree loaded.\n",
      "35: Tree loaded.\n",
      "36: Tree loaded.\n",
      "37: Tree loaded.\n",
      "38: Tree loaded.\n",
      "39: Tree loaded.\n",
      "40: Tree loaded.\n",
      "41: Tree loaded.\n",
      "42: Tree loaded.\n",
      "43: Tree loaded.\n",
      "44: Tree loaded.\n",
      "45: Tree loaded.\n",
      "46: Tree loaded.\n",
      "47: Tree loaded.\n",
      "48: Tree loaded.\n",
      "49: Tree loaded.\n",
      "50: Tree loaded.\n",
      "51: Tree loaded.\n",
      "52: Tree loaded.\n",
      "53: Tree loaded.\n",
      "54: Tree loaded.\n",
      "55: Tree loaded.\n",
      "56: Tree loaded.\n",
      "57: Tree loaded.\n",
      "58: Tree loaded.\n",
      "59: Tree loaded.\n",
      "60: Tree loaded.\n",
      "61: Tree loaded.\n",
      "62: Tree loaded.\n",
      "63: Tree loaded.\n",
      "64: Tree loaded.\n",
      "65: Tree loaded.\n",
      "66: Tree loaded.\n",
      "67: Tree loaded.\n",
      "68: Tree loaded.\n",
      "69: Tree loaded.\n",
      "70: Tree loaded.\n",
      "71: Tree loaded.\n",
      "72: Tree loaded.\n",
      "73: Tree loaded.\n",
      "74: Tree loaded.\n",
      "75: Tree loaded.\n",
      "76: Tree loaded.\n",
      "77: Tree loaded.\n",
      "78: Tree loaded.\n",
      "79: Tree loaded.\n",
      "80: Tree loaded.\n",
      "81: Tree loaded.\n",
      "82: Tree loaded.\n",
      "83: Tree loaded.\n",
      "84: Tree loaded.\n",
      "85: Tree loaded.\n",
      "86: Tree loaded.\n",
      "87: Tree loaded.\n",
      "88: Tree loaded.\n",
      "89: Tree loaded.\n",
      "90: Tree loaded.\n",
      "91: Tree loaded.\n",
      "92: Tree loaded.\n",
      "93: Tree loaded.\n",
      "94: Tree loaded.\n",
      "95: Tree loaded.\n",
      "96: Tree loaded.\n",
      "97: Tree loaded.\n",
      "98: Tree loaded.\n",
      "99: Tree loaded.\n",
      "100: Tree loaded.\n",
      "101: Tree loaded.\n",
      "102: Tree loaded.\n",
      "103: Tree loaded.\n",
      "104: Tree loaded.\n",
      "105: Tree loaded.\n",
      "106: Tree loaded.\n",
      "107: Tree loaded.\n",
      "108: Tree loaded.\n",
      "109: Tree loaded.\n",
      "110: Tree loaded.\n",
      "111: Tree loaded.\n",
      "112: Tree loaded.\n",
      "113: Tree loaded.\n",
      "114: Tree loaded.\n",
      "115: Tree loaded.\n",
      "116: Tree loaded.\n",
      "117: Tree loaded.\n",
      "118: Tree loaded.\n",
      "119: Tree loaded.\n",
      "120: Tree loaded.\n",
      "121: Tree loaded.\n",
      "122: Tree loaded.\n",
      "123: Tree loaded.\n",
      "124: Tree loaded.\n",
      "125: Tree loaded.\n",
      "126: Tree loaded.\n",
      "127: Tree loaded.\n",
      "Trees trained.\n",
      "Total training time:  469.8187255859375\n"
     ]
    }
   ],
   "source": [
    "# train random forest for Kaggle submission\n",
    "n_trees=128\n",
    "n_sample=len(train_X_max)\n",
    "n_feature=-1 # -1 ==> sqrt(d)\n",
    "max_depth=80\n",
    "\n",
    "rf = random_forest(n_trees=n_trees, n_sample=n_sample, n_feature=n_feature, max_depth=max_depth)\n",
    "\n",
    "start = time.time()\n",
    "rf.train(train_X_max, train_y) # use all data\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# predict and save labels\n",
    "test_y = rf.predict(test_X_max)\n",
    "df = pd.DataFrame({'Category': test_y.astype(int)})\n",
    "df.index.rename('Id', inplace=True)\n",
    "df.to_csv('spam_dist/spam_rf_tfidf_nt%d_sample%d_m%d_depth%d_nMax%d.csv' % (n_trees, n_sample, n_feature, max_depth, n_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. census"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(32724, 14)\n",
      "(16118, 14)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.preprocessing import LabelBinarizer, PolynomialFeatures\n",
    "from sklearn.feature_extraction import DictVectorizer\n",
    "train_data = pd.read_csv('./census_dist/train_data.csv').as_matrix()\n",
    "test_X = pd.read_csv('./census_dist/test_data.csv').as_matrix()\n",
    "train_y = train_data[:, -1]\n",
    "train_X = train_data[:, :-1]\n",
    "\n",
    "print(train_X.shape)\n",
    "print(test_X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# fill '?' with the most frequent element in that column\n",
    "def imputer(data, cols, method='max', missing='?'):\n",
    "    data_imputed = np.copy(data)\n",
    "    if method == 'max':\n",
    "        for j in cols:\n",
    "            col = data[:, j]\n",
    "            col = col[col != missing]\n",
    "            im = max_count(col)\n",
    "            data_imputed[data_imputed[:, j] == missing, j] = im\n",
    "            \n",
    "    if method == 'mean':\n",
    "        for j in cols:\n",
    "            col = data[:, j]\n",
    "            col = col[col != missing].astype(float)\n",
    "            im = np.mean(col)\n",
    "            data_imputed[data_imputed[:, j] == missing, j] = im\n",
    "            \n",
    "    if method == 'median':\n",
    "        for j in cols:\n",
    "            col = data[:, j]\n",
    "            col = col[col != missing].astype(float)\n",
    "            im = np.median(col)\n",
    "            data_imputed[data_imputed[:, j] == missing, j] = im\n",
    "    return data_imputed\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# impute data\n",
    "method = 'max'\n",
    "train_Xi = imputer(train_X, cols=np.arange(train_X.shape[1]), method=method, missing='?')\n",
    "test_Xi = imputer(test_X, cols=np.arange(train_X.shape[1]), method=method, missing='?')\n",
    "all_Xi = np.vstack((train_Xi, test_Xi))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# binarize categorical feature\n",
    "lb = LabelBinarizer()\n",
    "# all_Xib = \n",
    "# polynomials for the numerical features\n",
    "poly = PolynomialFeatures(2)\n",
    "all_Xib = poly.fit_transform(all_Xi[:, [0, 2, 4, 10, 11, 12]])\n",
    "\n",
    "for j in [1,3,5,6,7,8,9,13]:\n",
    "    b = lb.fit_transform(all_Xi[:, j])\n",
    "    all_Xib = np.hstack((all_Xib, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(32724, 126)\n",
      "(16118, 126)\n"
     ]
    }
   ],
   "source": [
    "train_Xib = all_Xib[:train_Xi.shape[0], :]\n",
    "test_Xib = all_Xib[-test_Xi.shape[0]:, :]\n",
    "print(train_Xib.shape)\n",
    "print(test_Xib.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(26179, 126)\n",
      "(6545, 126)\n"
     ]
    }
   ],
   "source": [
    "# divide data into train/validation\n",
    "S = np.random.choice(range(len(train_y)), int(len(train_y)/5*4), replace=False)\n",
    "train_data = train_Xib[S, :]\n",
    "train_labels = train_y[S]\n",
    "S = np.setdiff1d(range(len(train_y)), S)\n",
    "validation_data = train_Xib[S, :]\n",
    "validation_labels = train_y[S]\n",
    "print(train_data.shape)\n",
    "print(validation_data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total training time:  257.9419279098511\n",
      "Accuracy:  0.894457389511\n",
      "Accuracy:  0.838502673797\n"
     ]
    }
   ],
   "source": [
    "dt = decision_tree(15)\n",
    "start = time.time()\n",
    "dt.train(train_data, train_labels)\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)\n",
    "\n",
    "l = dt.accuracy(train_data, train_labels)\n",
    "print('Accuracy: ', l)\n",
    "\n",
    "l = dt.accuracy(validation_data, validation_labels)\n",
    "print('Accuracy: ', l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feature 54 <= 0.500000\n",
      "Feature 4 <= 7055.500000\n",
      "Feature 21 <= 488.000000\n",
      "Feature 12 > 1332.500000\n",
      "Feature 5 <= 2218.500000\n",
      "Feature 84 > 0.500000\n",
      "Feature 9 <= 394.000000\n",
      "Feature 1 > 32.500000\n",
      "Feature 8 > 11725338.500000\n",
      "Classification: 0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.])"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dt.predict(train_Xib[0, :], verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Trees = 20\n",
      "Trees trained.\n",
      "Total training time:  36.66999387741089\n",
      "Accuracy:  0.894457389511\n",
      "Accuracy:  0.852406417112\n"
     ]
    }
   ],
   "source": [
    "# train random forest\n",
    "n_trees=20\n",
    "n_sample=len(train_data)\n",
    "n_feature=-1 # -1 ==> sqrt(d)\n",
    "max_depth=80\n",
    "\n",
    "rf = random_forest(n_trees=n_trees, n_sample=n_sample, n_feature=n_feature, max_depth=max_depth)\n",
    "\n",
    "start = time.time()\n",
    "rf.train_parallel(train_data, train_labels, verbose=False) # use all data\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)\n",
    "\n",
    "l = dt.accuracy(train_data, train_labels)\n",
    "print('Accuracy: ', l)\n",
    "\n",
    "l = rf.accuracy(validation_data, validation_labels)\n",
    "print('Accuracy: ', l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n",
      "Feature 48 <= 0.500000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.])"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf.predict(train_Xib[0, :], verbose=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  8.96401572227478\n",
      "Accuracy:  0.76256684492\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  9.244059801101685\n",
      "Accuracy:  0.76256684492\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  8.947001457214355\n",
      "Accuracy:  0.805042016807\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  21.18677544593811\n",
      "Accuracy:  0.81436210848\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  24.4009530544281\n",
      "Accuracy:  0.812987012987\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  21.63522744178772\n",
      "Accuracy:  0.821390374332\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  41.40706014633179\n",
      "Accuracy:  0.838044308633\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  20.650776147842407\n",
      "Accuracy:  0.84782276547\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  21.616406440734863\n",
      "Accuracy:  0.852100840336\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  26.644749402999878\n",
      "Accuracy:  0.850267379679\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  42.28808116912842\n",
      "Accuracy:  0.850267379679\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  31.521824836730957\n",
      "Accuracy:  0.852864782277\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  25.180521488189697\n",
      "Accuracy:  0.853323147441\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  36.69718289375305\n",
      "Accuracy:  0.855614973262\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  31.006255865097046\n",
      "Accuracy:  0.853170359053\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  31.8244309425354\n",
      "Accuracy:  0.853628724217\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  35.67875671386719\n",
      "Accuracy:  0.854698242934\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  42.20683312416077\n",
      "Accuracy:  0.857142857143\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  33.989282846450806\n",
      "Accuracy:  0.856073338426\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  37.31991720199585\n",
      "Accuracy:  0.856073338426\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  41.34355902671814\n",
      "Accuracy:  0.858059587471\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  46.891753911972046\n",
      "Accuracy:  0.858823529412\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  36.329352378845215\n",
      "Accuracy:  0.853934300993\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  46.70114326477051\n",
      "Accuracy:  0.854545454545\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  43.0724093914032\n",
      "Accuracy:  0.857295645531\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  47.511680126190186\n",
      "Accuracy:  0.857601222307\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  44.98720932006836\n",
      "Accuracy:  0.858517952636\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  44.32628273963928\n",
      "Accuracy:  0.856684491979\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  50.88104200363159\n",
      "Accuracy:  0.854239877769\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  43.03969478607178\n",
      "Accuracy:  0.857906799083\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  45.190943002700806\n",
      "Accuracy:  0.859129106188\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  50.88291811943054\n",
      "Accuracy:  0.856073338426\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  50.88198494911194\n",
      "Accuracy:  0.856837280367\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  46.9473021030426\n",
      "Accuracy:  0.857295645531\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  65.7030930519104\n",
      "Accuracy:  0.854239877769\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  42.305405139923096\n",
      "Accuracy:  0.855309396486\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  49.113038063049316\n",
      "Accuracy:  0.857754010695\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  46.168540716171265\n",
      "Accuracy:  0.854851031322\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  53.577101945877075\n",
      "Accuracy:  0.857142857143\n",
      "# Trees = 50\n",
      "Trees trained.\n",
      "Total training time:  44.908695936203\n",
      "Accuracy:  0.856073338426\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff6eaf2de10>]"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# accuracy-depth relationship\n",
    "import matplotlib.pyplot as plt\n",
    "n_trees=50\n",
    "n_sample=len(train_data)\n",
    "n_feature=-1 # -1 ==> sqrt(d)\n",
    "\n",
    "all_depth = np.arange(1, 41)\n",
    "all_l = []\n",
    "for max_depth in all_depth:\n",
    "\n",
    "    rf = random_forest(n_trees=n_trees, n_sample=n_sample, n_feature=n_feature, max_depth=max_depth)\n",
    "\n",
    "    start = time.time()\n",
    "    rf.train_parallel(train_data, train_labels, verbose=False) # use all data\n",
    "    end = time.time()\n",
    "    print('Total training time: ', end-start)\n",
    "\n",
    "    l = rf.accuracy(validation_data, validation_labels)\n",
    "    print('Accuracy: ', l)\n",
    "    all_l.append(l)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAF5CAYAAABEPIrHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd4FNX+x/H3N40QIAlICUgVlWJDwN7FK8K1XkWIvYtd\nrNi9NqwgKCgWBCwIdhERQS82ihRB6agQehKEhJKElD2/P3bjL8QEUnYzm+Tzep48sLNnZr+HAfLJ\nOWdmzDmHiIiISLiI8LoAERERkaIUTkRERCSsKJyIiIhIWFE4ERERkbCicCIiIiJhReFEREREworC\niYiIiIQVhRMREREJKwonIiIiElYUTkRERCSshE04MbObzGyVmWWb2SwzO2Iv7S82swVmttPMNpjZ\nm2bWqFibBDMbHng/x8yWmdkZoe2JiIiIVEZYhBMz6wu8ADwCHA4sBKaYWeNS2h8HjAFeBzoDFwBH\nAq8VaRMNTANaA/8BDgSuBdaHrCMiIiJSaRYOD/4zs1nAbOfcbYHXBqwFhjnnni2h/Z1Af+fcAUW2\n3Qzc45xrHXjdH7gT6OicK6iCboiIiEgQeD5yEhjh6AZ8U7jN+RPTNOCYUnabCbQys16BYzQD+gCT\nirQ5K9BuhJltMrPfzOw+M/O8zyIiIlK6cPhG3RiIBFKLbU8FkkrawTk3A7gEGG9mucBGYCtwc5Fm\n++EPLBFAL+Ax/CMpDwSzeBEREQmuKK8LqAgz6wwMBR4FvgaaA88DI4FrAs0i8Aec6wIjMb+YWUvg\nLuDxUo67D9ATWA3khK4HIiIiNU4s0BaY4pz7qzIHCodwshkoAJoV294M2FTKPgOBn5xzgwOvF5nZ\njcAPZvaAcy4V/2hKrtt9Uc1SIMnMopxz+SUctyfwbkU7IiIiIlwMvFeZA3geTpxzeWY2D+gBfA5/\nL4jtAQwrZbc4ILfYNh/gAAu8/glILtamA7CxlGAC/hET3nnnHTp16lSOXlQ/AwYMYMiQIV6XEXLq\nZ82iftYsZe3n8rXpjPtxBsce2JHTu3WogsqCq7acz6VLl3LJJZdA4HtpZXgeTgIGA6MDIeVnYAD+\nADIawMwGAS2cc5cH2k8EXgtckTMFaAEMwX/FT+FoyyvATWY2DHgJ/6XE9wEv7qGOHIBOnTrRtWvX\n4PUuDCUkJNT4PoL6WdMEu5/bdu7ikxm/MvnXOTSIrcfIGy4jIsL2vmOI6XzC7KVreWbix3yz4UO2\nJf4E9R0T18TwR+5rvH7T5SXuE65qy/ksotLLIsIinDjnJgTuafIY/umcBUBP51x6oEkS0KpI+zFm\nVh+4Cf9akwz8V/sMLNJmnZn1xB9aFuK/v8kQ4B+XJotIzZeVk8fE2YuZ9Msc5myYS0reXLIb/AaR\neVAQBTvy2TUsl7G3X+t1qUHn8zkmzlrC+zO/Z3XGGh479yr+1e2Ave9YxaYv/JPnJn3Ed+kfsTNx\nNuTH0NT+RZ/Go7i19xlc+MpDvLH5Cubdt4AZ/32O2Jiw+BYmIRA2Z9Y5NwIYUcp7V5awbTgwfC/H\nnA0cG5QCRaRamTJ3BR/OnsnstXP4c9dcdtZfAFG7wBdBbP5BtIrqTteEq+l5SHfOO/ZQjnt8AG/n\n30LPb7tx8anV+6fc3LwCPpnxG+Nnf8fsTd+zMeZ7XN3NUBCF5Tfg9M+fpc17/RjW5wHOPrqzp7VO\nnrOcwV99yE9bPiI78RfIi6U5vbis2TsM/M+ZtG6a8HfbJc+8Rr/BXfjAdxv73rOIn+8ZT/sWjfZw\ndKmuwiaciIgEwxezl3L9+AfZkPAxOCMmtwP7RnSnS4Nk/nVQd84/tgtNG9b7x34zH36R5g/O5Yov\nL+DYTvNo17yhB9VXTE5uPuOmz+ejed8zJ+070mN/xMVmQH4dEnxHcVyd/px10Elc0eMYYqIjue7V\nN/k46xnO+epgWk64gCHnPcgFJxxapTU//cFUJi2ezsQvO0JuPVrZmVzQ8j7uOa8XSY3ql7hPRIQx\n4a6bGPxJZ+6a3YeOLxzJh+d/xjnHHlSltYeTqfNWctKh+xETHel1KUGlcFJLJScXXytcM6mfNcue\n+jlzyRouf+tRVsaNISqiNdc0Hs0jF55HyybxZTp2fL06TLn6A058tyvHPHs5657/lKhIb24FVdjP\nnNx8fl//Fys3pLMqLZ21f6WzPiOdtB3pbMnZTEZuOtt8m8isNx9idkBuHI3csZxS9w7O7XISl556\nJIn1Y/9x/Al33cSO7Gu5ceQYxu0YRJ9vDyPp43N47syHuKRHt5D37+F3vuDx5ecTd+D+3NbuKe46\n93Qaxdct8/53nHcK3dvPoefoczh30tHc9/s7PHXZOSGsuHJC8e9z+drNnPHi7ayOf5dGb5/Gj3e8\nR6fWTYL+OV4Ji9vXhwsz6wrMmzdvXm1bvCQSFG99/TPNEuM5o3uHKltYunRNOhe/8hS/RI3AchP5\nT+MHGXXjdcTXq1Oh4z367iT++/uZnBH1NJMfuDfI1e5ZVk4eZwx6nFk7xpMfk46L3frPRgXRROQ0\nJiavCXVdExpENqFzoy6c1/VELjq5G/XrxpT7M2978z3G/PkkefEraZLRiydOf4jrepV2g+7KuXf0\nJzz7Z19a7DiL5U+OK3e9RW3asoNuT1zOhoSP6WGP8/WDDwT1792atEx+Xp7CwtUpLNuUwuqMFDZm\nr2abbyPntr3ck/VJPp/jzlEfMHTlzRCRz8mxtzM962UifHUYedqHXN3zqKB91oc//EpeQQF9T+xS\npj/X+fPn061bN4Buzrn5lflshZMiFE5EKmbx6jR6DbuFtQkTAIja3pYOkb047+Be3HLmqSVOo1TW\nuvRtXDJ8MN/lvQAugh6xd/POzbeXOiVQHsc9/AAz7GmGHP4Nt597cuWLLYMfF62m1xsXsSPhZzpm\nX02b+HYkxTdh34ZNaNO4Ce2TmnDAvk1o2Tg+JMEvN6+AAW+O580VT7IrYQkNt/bgkVMe4rZzTgra\nZwx4YwIvrr2IVtsvYNlTbxMXG13pY+YX+Dj9iSf4H4/QMrMP8x56q1x/35av3cwHP81lzuolpGSm\nkJqzmgxS2BWb4p8aK1QQTdTO1jQoaEukRbM58Ss6Z/Vn9qNDKxWwymPBHxvp/fKNbEz8lBaZ/2Hy\nzcM5dL8k5q5Yzykj+rAjfi59E17kvQE3VOrvyNI16Zw97B5+bzAagIidSeznO4MzO/bi1n//q9Qp\nT4WTEFE4ESkfn89x82vv8erq2wDj+jZDaBbfkI9/m8zS/MnkN/gT8mNotP1Ejk/qxbWn9KL3ER0r\n9R9nxo4crnj5FT7PeBIXtZPu7mbeu2EgB7TcJ2j9ysnNp/k9p7OtzhLmXfcLXdo3D9qxS3L3qI94\nfuXVROYlMqLHuJCNWpRFfoGPe0d/zIhFj5OT+Cv7ZPTk06tf5fiD21bquDe++i6vbLyMtjuSWTpo\ndNCvtLl39Cc8u/JSYrP2Z+rVn5ZYb0pqBhN+nMe3y+by219z2RQxh4IGKf438+oSk9WWeNeGJtFt\naB3flgOatOHgVm044oC2HLpf0m7TfJe9+Dpvb7mJhG3HMPOOD0M6peLzOa4ZPprRG+7AfHW4o8Nw\nnrvq/N3a7MjO5fjH7mZh7DDabruY2Q+OLPcPBfkFPq546Q3eS/Vf+JrcdBBd2xzA+/Mm81v2ZHYl\nLAZfJA0yj+aofXpxxXG96HtSl7//XBROQkThRKTs5ixfx5mv9ictcRKtt/Xjq9uG7fYftM/nmDp/\nJSO/ncz3GyfzV/3pELWLyO1t6BDZi3MP7sX1PU+iUYOyrTXIzSvgztHjeHvtoxTEbaBj9lWMvfph\njujQMiT9W7QqlcNeOZwGuQey6dlpIblsdcu2bE548g6WxL3KvpkX8NO9r9OmWWLQP6ci8gt8PPD2\nZ7yw+DYKYrZwfuIg3r/jpgqtw7l2+BjeSL+S/XdczuJBb4Rs8eZHP/5Gv0/PoSBqO48f/g7xdesy\nbclcFqTPYSNzyYv/3d8wtz4JO7uxf73uHNumO+cc0Z1TDmtf7tA84osfufmH84nwxfLumZ/R96Qu\nQe/Tj4tWc84b17Gl4VT2234ZX98xZI9XKN362vu8lHINdbLb8tlFH9Oz+4Fl+pxx03/h2s9uYGfi\nbPbffgWf3/rsPwLXzCVrGP71V3ybMpmNcdMgZgcRWc3Yr6AXZ3bsxamtmnL2GaeAwklwKZyI7F3h\nT1fvpt1NRH59Bh7yCk9eevZe99ucmcXwSdP56NfJLM2bTH78HxX6/FaZF/L6RY+X+T/dynh54g/c\nMvcUjiq4i1lPPB3UY38+awkXTujLrrjfubjxUMbeem1Y3ACuuHXp2zjj+YEsjnuFBluPY/zFb9Lr\niLLfpfXyoW8wdut1dMy6ht8GvRryRcYr1/3FUc9fyNaG3/o35NWlwY6utK/bnaNbd+esbt05vduB\nQatj5pI19Hj9PLLjljFgv9EMvrpPUI6bX+AjefAIPswYSGRuIx7uOpKHk3uVad/PZy2hz4TzyY1d\nz137v/WPUZai1qVv46zBD7Mg5iXqbO/E86e8ws1nnbDXz9iRnctrX/3E+3Mn82v2ZHYlLIL1Bq87\nUDgJLoUTCQfPfjiNKUtm8u7NA4KyfiKYvl3wB+ePvoaMhtPpsOMaptz1XIV/0p86byWfzPmZ/IKC\nMu9z2iGHceGJh1Xo8yrqzEHPMyn3bu7f77MyhbC98fkcV770JmPTb6VOVjvGnT+e8447OAiVhtbQ\nz77j7u+vIS9uLT3rPMqn99y119Gk5Bde4f0dN3JI9k3Mf3JYlV39lJWTx3MfT+Xg1q3495GdQn6z\nts2ZWXR77GrWxL/P8b4H+d/D/61UXyfPWU7fd69me8OfOCjrBr666+kyX3VWaMNf2znqyatZl/AB\n3XLv5MdHnt7tz8Hncwx4cwIvrxyAL3ob/67/KBMG3FbhdUCzl67lsTdG8uXgJyEI4QTnnL4CX0BX\nwM2bN8+JVLV5K9a7VgP6Oh7F8Yi5yDvbuKfGT/G6LOecc7ty893Zg15wPFDXRd3Zzj374TSvS6oy\nBQU+1/z2c50NTHT/W/BHpY6Vkprx9znuePe1Lj1jZ5CqrBrpGTtd9/vvcjwc4ere3tW9P/2XUtue\n98yLjkdxhw+83RUU+KqwSm8UFPjcGU887XjEXNLtZ7u1aZnl2n971i73/EffuMMH3u54sI6LvmN/\nN+yz7ypd07lPD3E8FOUSbjvR/fL7Buecc1/NWe4a3Xaa41FciwH/cbOWrKnU5xSaN2+ew/+Mu66u\nkt+PNXJShEZOxAs5uflc9OJwPsl4CCuI5bo2L9DvuOP4z1vXsbXhN+y//Qq+vnOwZzcF+2zGYi7+\n4Cp2Jsyhy65bmXLPkyG5+iacpaRmsP+z3Yn2JbDh8Z9KvHfI3rz19c9c93U/8mP+4vb9XmfINReG\noNKqMWbqHK6ffBW76i/jOO7ly3sf2u3S7cLRpiPz7mHmY0+H5XRVqDz67iT+u/gi6uS05MvLPufU\nLu1LbTtzyRpenjKZb9dMZlPcN/41HDuTOKrOlXxx90PluvfLnoz44kdu+f5CwHFI5IUsjHqVqOx9\nebDryzxyUe+gfAZoQWzIKJxIVXvjq1ncOuUGshMWclB2fybe/uTfIeTv4f/UO4nIj+OuziN45orz\nqqSuNWmZDJ04lYnLvmRl3XeI2dGeYae9yfW9a+/TIMZ/t4B+046m467LWfrsyDLtk5tXwI+LV/PM\nFxP4Ovdh6m3vypdXvs+Jh7YLcbWhtyM7l7OeHcT0gieJ2bE/w08fxTVnHM3pjz/FVN8DHO97kO8e\neaxWBZNCX/68jHPHn0N+TDpPd5vAPRecBvgfNDnyqx8ZP38yi7InsythCfgiic88lqMb9+KqE3rR\n54TDQvJn9uufmzhpWDIZDWZwQsS9fHrnfUELP4UUTkJE4USqyh8btnDmi/exLO514jK7MqL3K1z+\nryNKbDt3xXrOfOUGUhMn0jKzD5NveYmD2zULaj0+n+ODHxYy6ofJzNo8mW0JMyCigJjMThyfmMxH\nd9xdodGCmuaKYW8yZus1XNt4DK/ddNnf27Ny8pj+6x9MX7yE+WuX8HvmUtJ8S8iutwyi/Q9oPTLv\nHv734BNBubdHOPnkp0Vc+uFV7EyYS8OMU9na8BtOtcf45uGHvC7NUympGXR/OpnN8V9zeO6tbMz5\nMzA6spOInc1p787grI69ufXM06rsCq38Ah8pqRkhex6RwkmIKJxIqOUX+LhuxBhGr78HF5FLn0ZP\n8c5t/fd6aaXP57j9jfG8/OctgI/r2wxl+PUXV+onrJTUDIZ+MZUvlk3mD/sKX72NkFuPpKwenNKq\nFzf17MVxB7Wp8PFrqgPvvoqVdd7naLuNtTt/ZzNL2FV/pf/pxoDlNKRB9kG0iOlMx306cWS7zvQ8\n/GC6HtDC48pDJyc3nwteeJFJO/5Lz7oP8NWDA/e+Uy2Qm1fA8Y/exxx7ifgdR4R8dMRrCichonAi\nofTRj79x1Yc3sq3hj7TddjGf3fQ8h+6XVK5jLF2TTq+ht5ESP44mGb2ZeN2rHNWp1R738fkcS9ek\nM3dlCgvXrObX9cuZs2UK2xJmQkQBdTI7c3DdXvTt2ovrzzi+wrd9ry02Z2ax/yO92R6znPjcTrSK\n7UznJp05un1nTjusE53bNK2R33jKIr/A59nziMKZz+dqxd+JYIYTPfhPJARycvOZu2Id8/5IYfG6\nFGat/Znf6rxKTMQBvHDot9xx3ikVOm6n1k1Y/cJ7PPR2Mk/91p+j3z6Ii5o8x+29/83c31NYtDaF\n3zensGbbatLzUthuKeTGpUB09v8fZFc8SXYKZyYO5+aevTimc+sg9bp2aJwQR8aL070uIywpmJSs\nNgSTYFM4EamAnNx8fvhtFb+sSmHx+tX8uSWFDTtT2Jy/mp3RKRTErYeI/79/R4Rrxr+iHuPjJ+8K\nynM4Hr/0LK5JPYGez9/Ne9v68977//+eZTciNq8NibShc52etG3QlgObteHQ1m048sA2tG/RSP9Z\nikhYUzgRKacxU+dw3VeXkhu/3L/BGRFZLaiX34ZGkW05KPZ42iW2oXOLthzerg1HdmhN44S4oNfR\nplkiy557nTFTr2NVWhqHtW3DUR3b0GKfBkH/LBGRqqRwIlJGWTl59H76Sb5zTxDnO5yHO37FMR32\n58gOrarsqaQlKe0qHxGR6krhRKqNO978gG9W/sDQ5Ns5+bD9qvSzJ89ZTp/3LmVn/HxOtAeYPOjB\nGndJqIhIuNDqJakWBk34miEpF/Fr1Ouc8vGB7H/XFUyZuyLkn5tf4KPPcy/T+7PDyY3IZNRxM/ju\n0f8qmIiIhJDCiYS9Cd8v5P4FF9Bk++msH5DGefVeYFXEVM6Y2Im2d17EZzMWh+Rz565YT7M7z+DD\nrFs4JP8qNjz8C1eefmRIPktERP6fwomEtdlL15L8RW/qZh/Arw+Np8U+Dfj4ntv46+E/6Bf/Musi\nfuLcqQfT8o4LGP/dgqB97q2vvc+Rbx1CRswinuj4Fb8+/XJIFrWKiMg/KZxI2FqTlslJr/fGfNHM\numUSSY3q//1eYv1Yxt15AxmPreSyxNdJjfiFftMPJ2nA2YyZOqfCn/nHhi20uTOZlzYm03LX6ay4\nfREP9O0ZjO6IiEgZKZxIWNqRnUuXp84nt846Pu0zudQ7qdavG8OY265h+1PLubbxGLZELOeKGUfS\nZEAvRn45gx3ZuWX+GjThaw4ccihrY77i5qT3WDP4/ZA9g0JEREqnq3Uk7Ph8jsMeuoatDX5g2BFT\nOfOoTnvdJzYmitduuoyX8y7mrrc+4LXMx+k/5zj6l3MQpVHeaXx1zVsc0aFlBasXEZHKUjiRsHPS\nfx/mzwZvc0vzcdxy9onl2jcmOpJh1/VjcMGFPPPhVFI2p5Z536SERjyc3Fu34BYR8ZjCiYSVy4e+\nwY8RT9Ar6hmGXdevwseJiozQWhERkWpK4UTCxhPvf8XYLf05KOcGvhh0t9fliIiIRxROJCyMm/4L\nD/3Wh6ZZvZj7zDA9mE5EpBZTOBHP/bQ4hUu+7E1cXkd+e/h9YmP011JEpDbTyj/x1KqNW+nxZm8i\nfHWZfesXNG1Yz+uSRETEY/oRVTyTlZNH12f+Q27sJiadP4OD2zXzuiQREQkDCifimYFjPyKj4XRe\n7DKdXkd08LocEREJEwon4pl3l4yiAcdx2zkneV2KiIiEEYUT8cRPi1PYkjiNK/d50+tSREQkzGhB\nrHjioQ/HQF4cT13cx+tSREQkzGjkRKpcfoGPH7a/xQHWd7cnDYuIiIDCiXhgyKf/I7/BagZ0v8rr\nUkREJAwpnEiVGz5jFNEcyPW9jvW6FBERCUNacyJVatXGraTEfUSPRlfpFvUiIlIihROpUve/9z5E\n5DOo32VelyIiImFK0zpSpT5fN4pm9KZL++ZelyIiImFK4USqzIc//EpW4lxubfuJ16WIiEgYUziR\nKvPkl6MwX1Me6vtvr0sREZEwpjUnUiW27dzFQvcOXaMuJS422utyREQkjGnkRKrEf9+fiKv7Fw/3\n1L1NRERkzxROpEq8vWgU9Tmas4/u7HUpIiIS5jStIyE3Z/k60uOncF4bjZqIiMjeKZxIyN0/YQwU\n1OGpi/t6XYqIiFQDmtaRkMov8DE9cxT70YeWTeK9LkdERKoBhRMJqeFf/EB+gz+5tcsor0sREZFq\nImymdczsJjNbZWbZZjbLzI7YS/uLzWyBme00sw1m9qaZNSqlbT8z85nZx6GpXkoz7IdRRG1rzy1n\nneh1KSIiUk2ERTgxs77AC8AjwOHAQmCKmTUupf1xwBjgdaAzcAFwJPBaCW3bAs8B34egdNmDdenb\n+DP2A05J1EP+RESk7MIinAADgJHOubHOuWVAfyALKO3yjqOBVc654c65FOfcDGAk/oDyNzOLAN4B\nHgZWhax6KdH9746HyF08eaEe8iciImXneTgxs2igG/BN4TbnnAOmAceUsttMoJWZ9QocoxnQB5hU\nrN0jQKpz7q1g1y1793HKmzTZ1pMjOrT0uhQREalGwmFBbGMgEkgttj0V6FDSDs65GWZ2CTDezGLx\n9+Nz4ObCNmZ2PHAlcFgoipY9+2zGYnYmzub6lh94XYqIiFQzno+cVISZdQaGAo8CXYGeQDv8UzuY\nWX1gLHCtc26rR2XWao9PegvL3odH+p3ldSkiIlLNhMPIyWagAGhWbHszYFMp+wwEfnLODQ68XmRm\nNwI/mNkDQBLQBphoZoUrMSMAzCwX6OCcK3UNyoABA0hISNhtW3JyMsnJyWXvVS2WlZPH/PyxHBZx\nCfH16nhdjoiIBNm4ceMYN27cbtsyMzODdnzPw4lzLs/M5gE98E/NEAgUPYBhpewWB+QW2+YDHGDA\nMuCQYu8/CdQHbgXW7qmmIUOG0LVr13L0Qop6fPwkXFw6D/a42utSREQkBEr6gX3+/Pl069YtKMf3\nPJwEDAZGB0LKz/iv3okDRgOY2SCghXPu8kD7icBrZtYfmAK0AIYAs51zhaMtS4p+gJll4F9ruzTE\nfan13lowiji6c/7xxfOhiIjI3oVFOHHOTQjc0+Qx/NM5C4Cezrn0QJMkoFWR9mMC60puAp4HMvBf\n7TOwSguXf1jwx0ZS47+kX/xLXpciIiLVVFiEEwDn3AhgRCnvXVnCtuHA8HIc/x/HkOC77/2xUBDN\n05dofY6IiFRM2IQTqRmmb36P1hHn0qZZoteliIhINVUtLyWW8LQmLZOchN/4V7ueXpciIiLVmMKJ\nBM2702eDOc4/qrQb+4qIiOydpnUkaKYsmYHlNKJntwO9LkVERKoxhRMJmkUZM2lix+oJxCIiUima\n1pGgyM0r4K+4WRzWSFM6IiJSOQonEhQTZy+BOtvofcixXpciIiLVnMKJBMUnc2eCL5KLTjrC61JE\nRKSa05oTCYpZ62ZQt+Awmjas53UpIiJSzWnkRIJijZvB/nU0pSMiIpWncCKVtnztZvLiV3JSO4UT\nERGpPIUTqbR3ps8EoO9xulJHREQqT2tOpNKmrZhBxK7mHNu5jdeliIhIDaBwIpW2dPtMmkfo5msi\nIhIcmtaRSsnKySOz/s8c3kRTOiIiEhwKJ1IpH/64EKKzOauLFsOKiEhwKJxIpXz+y0zIj6HfiV29\nLkVERGoIrTmRSpmzaQb1C7oRX6+O16WIiEgNoZETqZT1ETPoUE9TOiIiEjwaOZEKm7tiPQX113Bq\na4UTEREJHo2cSIWN+8F/87WLTtSVOiIiEjwaOZEKm/7HDKJy29KlfXOvSxERkRpEIydSYSuyZ7Kv\n05SOiIgEl8KJVEjGjhx2NJjHEUma0hERkeBSOJEKGffdPIjM49xuGjkREZHgUjiRCpn060zIjeP8\n4w/1uhQREalhtCBWKmR+2gwS3ZHExuivkIiIBJdGTqTcfD5HaswMOsdrSkdERIJPP/ZKuf24aDW+\nuFR6tNdiWBERCT6NnEi5vT9jBgAXn3S0x5WIiEhNpJETKbcfVs0gJr8DHVo19roUERGpgTRyIuX2\nR+5MWpumdEREJDQUTqRcNm3ZQXb8Qo5pqcWwIiISGgonUi7vTv8ZInycf6TCiYiIhIbCiZTLV4tn\nQk4C/z6qk9eliIhIDaUFsVIuC7fMYB+OJipSuVZEREJD32GkzPILfGyOncmhiZrSERGR0FE4kTL7\net4KXOxWTu+sK3VERCR0FE6kzD6YNQOcccnJR3ldioiI1GBacyJlNmPtDGLzD6Flk3ivSxERkRpM\nIydSZikFM2kXrSkdEREJLYUTKZNVG7eyK2EJx7fRYlgREQkthRMpk7enzwKg77EKJyIiElpacyJl\nMm3ZTGxXY045rL3XpYiISA2ncCJlsihzBs3sWCIizOtSRESkhtO0juxVbl4BW+vNpktjTemIiEjo\nKZzIXn06cxHE7KD3IbpSR0REQk/hRPbqkzkzoCCK5JO6e12KiIjUAlpzIns1e8MM6vkOp3FCnNel\niIhILaCRE9mrtcxk/1hN6YiISNXQyIns0aJVqeTH/8FJLbQYVkREqkbYjJyY2U1mtsrMss1slpkd\nsZf2F5vdzat7AAAgAElEQVTZAjPbaWYbzOxNM2tU5P1rzOx7M9sS+Jq6t2PKP737/UwALjpe4URE\nRKpGWIQTM+sLvAA8AhwOLASmmFnjUtofB4wBXgc6AxcARwKvFWl2EvAecDJwNLAW+NrMmoemFzXT\ne7+OJ3rbARzVqZXXpYiISC0RFuEEGACMdM6Ndc4tA/oDWcBVpbQ/GljlnBvunEtxzs0ARuIPKAA4\n5y51zr3qnPvVObcCuAZ/f3uEtCc1yMwla1hT/wPOSbrV61JERKQW8TycmFk00A34pnCbc84B04DS\nVmHOBFqZWa/AMZoBfYBJe/ioekA0sCUIZdcKt7w7DMttwEtXX+F1KSIiUot4Hk6AxkAkkFpseyqQ\nVNIOgZGSS4DxZpYLbAS2Ajfv4XOeAdbjDz2yF+vStzHPvc6RkdeT1Ki+1+WIiEgtEg7hpNzMrDMw\nFHgU6Ar0BNrhn9opqf1A4ELgXOdcbhWVWa3dMmoURGXx0qV7ynsiIiLBFw6XEm8GCoBmxbY3AzaV\nss9A4Cfn3ODA60VmdiPwg5k94Jz7exTGzO4C7gF6OOcWl6WgAQMGkJCQsNu25ORkkpOTy7J7tZeT\nm88X6UNp5/pxRIeWXpcjIiJhZty4cYwbN263bZmZmUE7frnDiZnt55z7M1gFOOfyzGwe/oWqnwc+\nwwKvh5WyWxxQfATEBzjg78fmmtk9wH3A6c65X8pa05AhQ+jatWuZ+1DT3P/2J+Q3WM3jJ3zkdSki\nIhKGSvqBff78+XTr1i0ox6/ItM7vZvY/M7vEzGKDUgUMBq41s8vMrCPwKv4AMhrAzAaZ2Zgi7ScC\n55tZfzNrF7i0eCgw2zm3KbDPvcBj+K/4WWNmzQJf9YJUc431+qLBJG49mYtPrb0BTUREvFORcNIV\n+BV/oNhkZiPN7Mi97LNHzrkJwF34w8QvwKFAT+dceqBJEtCqSPsxwB3ATcBvwHhgKXB+kcP2x391\nzofAhiJfd1am1prutckz2ZE4i5u73eF1KSIiUkuZ/6rdCuxoFgWcDVwBnAGsAEYBbxcJFdWKmXUF\n5s2bN6/WTuu0vKMPafYrWc8uJSqyWq6XFhERDxSZ1unmnJtfmWNV+LuPcy7fOfcx/vuL3AvsDzwP\nrDWzsboTa/Xz/a+rWN/gYy7Yd4CCiYiIeKbC34HMrLuZjcB/j5E78AeT9sC/gBbAZ0GpUKrMbeOG\nYrsSGXb1ZV6XIiIitVhFrta5A7gS6AB8CVwGfOmc8wWarDKzK4DVQapRqkBKagYLIt7kuMjbaJwQ\n53U5IiJSi1XkPic34F9bMto5t7GUNmnA1RWuSqrcTaPegIhcRlyhm66JiIi3yh1OnHMHlKFNLv6n\nBks1kJWTx1dbhrK/u4hD9yvxiQEiIiJVptxrTszsSjPrU8L2PmZ2eXDKkqp075gPKai/jkHnDPC6\nFBERkQotiL2Pfz6kD/xTOfdXrhypaj6fY9SywTTaehoXnHCo1+WIiIhUaM1Ja2BNCdtTAu9JNTJi\n0o9kJc7lvg6TvS5FREQEqNjISRr+O7gWdxjwV+XKkar25LcvUCezM/df2NPrUkRERICKjZyMA4aZ\n2Xbg+8C2k/A/2+b9YBUmoTd13ko2JXzOZQ1fIyLC9r6DiIhIFahIOHkIaAt8A+QHtkUAY9Gak2rl\nzg+GYq4xQ++6xOtSRERE/lbuaR3nXK5zri/QEbgY+A/Q3jl3VeASYqli17w8msMG3sqbU2bj85Xt\nWUl/bNjCb5FvcVLcTSTWD9bDpUVERCqvIiMnADjnVuB/2J94aMu2bEatvwMXvZNrZr3EDV8fyEmJ\nl/Hofy7huIPalLrfjaNGQkQBL19+QxVWKyIisncVCidm1hL/E4lbAzFF33PO3RGEuqSM7n/nQ1zs\nVr7svYxFa9bx6syxTMsZxLQPHyTh9ZM4b7/LeDz5Alo2if97nx3ZuUzb9hId3KUc1Laph9WLiIj8\nU0VuwtYDWI7/NvZ3Aqfgf9bOVUCXoFYne/XeipE03NqDXkd04O7ze/DH82PYeOcmrm8ylkiiGb3l\nGloNbUbbOy/iife/Iic3n7veGo+v3kaeOe92r8sXERH5h4pcSjwIeN45dwiQA5wPtAK+Az4IYm2y\nF5/8tIjtDX/iikOu3217UqP6vHrjpfz14lR+Tl7DGXUfZZNbyEPLe1HvwVa8vuo+GmecwTnHHuRR\n5SIiIqWrSDjphP/KHPBfrVPXObcDeBi4N1iFyd49MnEkltWUxy46p9Q2R3RoyeQH7iXr+UW8ffxc\nDo3sS6SvLo+d9mAVVioiIlJ2FVlzspP/X2eyEWgPLA68bhyMomTvNmdm8Zu9zdHRN1C/bsxe20dE\nGJf06MYlPboBL4a+QBERkQqqSDiZBRwPLAW+BF4ws0PwX1I8K4i1yR7cM3Y8xGYy6NxrvS5FREQk\nqCoSTu4A6gd+/0jg932BlYH3pApM+HMk+7jTOfmw/bwuRUREJKjKFU7MLBJoCfwK4JzbCfQPQV2y\nBxO+X8jOxNnc2Pojr0sREREJunItiHXOFQBfAw1DU46UxWOTRhKxM4lHk8/yuhQREZGgq8jVOosA\nzSV4ZNOWHSyOfIdjYq8mLjba63JERESCriLh5EHgeTM708yam1l80a9gFyi7u/ft9yFmB8/01UJY\nERGpmSqyIPbLwK+fA0WfMmeB15GVLUpK9+HqkTThjD0+N0dERKQ6q0g4OSXoVUiZvPvtfLIS53Jb\nu0+9LkVERCRkyh1OnHPfhaIQ2bsnvhpJhNuXh/v92+tSREREQqbc4cTMTtzT+8657ytejpRmw1/b\nWRb9HidG30FsTIUeJi0iIlItVOS73PQSthVde6I1JyFw99j3ICqLZ/te43UpIiIiIVWRq3UaFvtq\nCpwBzAFOD15pUsjnc3yyZiRNt/XmqE6tvC5HREQkpCqy5iSzhM1TzSwXGAx0q3RVspu3v5lLduIv\n3NX+Ma9LERERCblgLl5IBToE8XgS8NTXI4mkFQ/27eV1KSIiIiFXkQWxhxbfBDQHBgILglGU/L81\naZmsiBnHKTH3EhOt5TwiIlLzVWTkZAH+BbBWbPss4KpKVyS7ueftdyFqF88mX+11KSIiIlWiIuGk\nXbHXPiDdOZcThHqkCJ/P8dm6kSRxJt0P3NfrckRERKpERRbEpoSiEPmnUV/PJifxV2448GmvSxER\nEaky5b6U2MyGmdnNJWy/2cxeDE5ZAvDMtJFEbm/DwAt0hbaIiNQeFbnPyfnAjyVsnwFcULlypFBK\naga/x47nlIRrtRBWRERqlYqEk32A7SVs3wY0rlw5UujOsW9DRB7PXaQ1xiIiUrtUJJz8DpR0w41e\nwJ+VK0fAvxD2i40jabHjbLq0b+51OSIiIlWqIlfrDAZeNrMmwLeBbT2AO4Hbg1VYbfbW1J/ZlbCY\nmzsN9roUERGRKleRq3VGmVkd4AHgocDm1cANzrmxQayt1vppxRIAbjnrZG8LERER8UCFbl/vnHsF\neCUwepLtnNsR3LJqtw2ZqVhOQ+rXjfG6FBERkSpXkdvXtwOinHMrnXPpRbYfAOQ551YHsb5aKW1n\nKtG7mnldhoiIiCcqsiB2NHBUCduPCrwnlbQlN426PoUTERGpnSoSTg4HZpawfRbQpXLlCMC2glQa\nRDT1ugwRERFPVCScOCC+hO0JgO4WFgTZEWk0jNHIiYiI1E4VCSffA/eZ2d9BJPD7+yj5zrFSTrui\nU2kSp5ETERGpnSpytc69+APKcjP7IbDtBPwjJ6cEq7DaKjevABe7mebxGjkREZHaqdwjJ865JcCh\nwASgKdAAGAscGNzSaqeV6/+CCB+tGymciIhI7VTR+5xsAO4HMLN4oB/wFdAdrTuplGXrUgFo20TT\nOiIiUjtVZM0JAGZ2opmNATYAdwH/A46uxPFuMrNVZpZtZrPM7Ii9tL/YzBaY2U4z22Bmb5pZo2Jt\n+pjZ0sAxF5pZSc8ECit/pqYBcGALjZyIiEjtVK5wYmZJZjbQzFYCH+B/EnEd4Fzn3EDn3JyKFGFm\nfYEXgEfwX6q8EJhiZiU+5djMjgPGAK8DnYELgCOB14q0ORZ4L9CmC/AZ8KmZda5IjVVldbp/5KRj\nK42ciIhI7VTmcGJmE4Hl+Neb3A60cM7dEqQ6BgAjnXNjnXPLgP5AFnBVKe2PBlY554Y751KcczOA\nkfgDSqFbgcnOucHOueXOuYeB+cDNQao5JNZnpEJuHEmN6ntdioiIiCfKM3LSC3gTeMQ5N8k5VxCM\nAswsGugGfFO4zTnngGnAMaXsNhNoVThNY2bNgD7ApCJtjgkco6gpezhmWEjdmUaUbl0vIiK1WHnC\nyfH4r8yZZ2azzezm0qZdyqkx/kW0qcW2pwJJJe0QGCm5BBhvZrnARmAru4+KJJXnmOHir5xU6hRo\nSkdERGqvMocT59ws59y1QHP8Uyj98C+GjQD+ZWYNQlPiPwXWjQwFHgW6Aj2BdoG6qrXM/DQaoJET\nERGpvcp9KbFzbicwChhlZh2Aq4GBwNNmNtU5d3Y5D7kZKIB/fEduBmwqZZ+BwE/OucGB14vM7Ebg\nBzN7wDmXGti3PMf824ABA0hISNhtW3JyMsnJyXvbtdJ2kEqrKD2iSEREwte4ceMYN27cbtsyMzOD\ndvwK3eekkHNuOXCPmd0HnEXpC1j3dIw8M5sH9AA+BzAzC7weVspucUBusW0+/M/9scDrmSUc41+U\n/NDC3QwZMoSuXbuWtQtBtSsqlcZ1NXIiIiLhq6Qf2OfPn0+3bt2CcvxKhZNCgcWxnwa+KmIwMDoQ\nUn7Gf/VOHDAawMwG4b866PJA+4nAa2bWH/8i1xbAEGC2c65wZGQoMN3M7sC/UDYZ/8LbaytYY8j5\nfI6C2DSSGiiciIhI7RWUcFJZzrkJgcW1j+GfelkA9HTOpQeaJAGtirQfY2b1gZuA54EM/Ff7DCzS\nZqaZXQQ8GfhaCZwTuP1+WFq3eRtE7WLfRC2IFRGR2isswgmAc24EMKKU964sYdtwYPhejvkR8FFQ\nCqwCS9f4Ly5q10QjJyIiUntV+Pb1Enx/bPLfun6/JI2ciIhI7aVwEkZWpflHTjq11MiJiIjUXgon\nYWTd1jQoiKJd84ZelyIiIuIZhZMwsnF7KhE5TYiK1GkREZHaS98Fw0h6dip18jSlIyIitZvCSRjJ\nyEsjzmkxrIiI1G4KJ2Fkuy+VhCiNnIiISO2mcBJGciLTaFRH4URERGo3hZMwklcnlab1NK0jIiK1\nm8JJmMjYkQN1trFvgkZORESkdlM4CRNL1/jvDtt6H42ciIhI7aZwEiZWrPffHbZ9kkZORESkdlM4\nCROr0vwjJx32VTgREZHaTeEkTKT85R856diqiceViIiIeEvhJExs2JaKZTciLjba61JEREQ8pXAS\nJtJ3phGdq8WwIiIiCidhYktuKnV9Wm8iIiKicBImthWkER+hcCIiIqJwEiayIlJJjNG0joiIiMJJ\nmMiNTqVpnEZOREREFE7CQG5eAa7uZpLiNXIiIiKicBIGVqzbDOZo3UgjJyIiIgonYWD5ev/dYfdr\nqnAiIiKicBIG/tjkvzvs/s01rSMiIqJwEgZWb/aHk86tNXIiIiKicBIGNmSkQW4cTRvW87oUERER\nzymchIFNO1KJ2qVRExEREVA4CQtbdqURW6BwIiIiAgonYSEjP5X6aDGsiIgIKJyEhZ2kkhitkRMR\nERFQOAkLu6LSaFxX4URERAQUTjzn8zkKYlNpVl/TOiIiIqBw4rl1m7dBVC6tGmrkREREBBROPLd0\njf8GbG0aa+REREQEFE48t3Jj4a3rNXIiIiICCieeS0n3P/SvUyuFExEREVA48dzaralQEEWbZole\nlyIiIhIWFE48tml7GhE5TYmK1KkQEREBhRPPpWenUidPi2FFREQKKZx4bGtuKnFovYmIiEghhROP\n7XBpJEQqnIiIiBRSOPFYTmQqjepoWkdERKSQwonH8uqk0ay+Rk5EREQKKZx4KGNHDtTZRot4jZyI\niIgUUjjx0OKUwlvXa+RERESkkMKJh37f4L877P5JCiciIiKFFE489Eeqf+TkgBaa1hERESmkcOKh\ntVv8IycdWzXxuBIREZHwoXDioQ3bUrHsRsTFRntdioiISNhQOPFQ2s5UonO13kRERKQohRMPbc1N\nI86ncCIiIlJU2IQTM7vJzFaZWbaZzTKzI/bQ9i0z85lZQeDXwq/firW73cyWmVmWma0xs8FmVif0\nvSmbbQWpNIjQYlgREZGiwiKcmFlf4AXgEeBwYCEwxcwal7LLrUAS0Dzwa0tgCzChyDEvAgYFjtkR\nuAq4EHgyNL0ov6yIVBrGaORERESkqLAIJ8AAYKRzbqxzbhnQH8jCHyj+wTm33TmXVvgFHAkkAqOL\nNDsG+NE5N945t8Y5Nw14P9A2LORGp9EkTiMnIiIiRXkeTswsGugGfFO4zTnngGn4A0ZZXAVMc86t\nLbJtBtCtcHrIzPYDegOTglF3ZeXk5uPqbqZFvEZOREREioryugCgMRAJpBbbngp02NvOZtYc6AX0\nK7rdOTcuMC30o5lZ4DNedc49E5SqK+n39X+BOVrto3AiIiJSlOcjJ0FwBbAV+KzoRjM7Gbgf/xTR\n4cB/gDPN7MEqrq9Ey9b5s1i7JprWERERKSocRk42AwVA8SGEZsCmMux/JTDWOZdfbPtjwNvOubcC\nrxebWX1gJPDEng44YMAAEhISdtuWnJxMcnJyGcopm983+cPJgS00ciIiItXLuHHjGDdu3G7bMjMz\ng3Z8z8OJcy7PzOYBPYDPAQLTMD2AYXvaNzA60h54s4S344DigcVXePzAupYSDRkyhK5du5a1CxWy\n5q/CW9dr5ERERKqXkn5gnz9/Pt26dQvK8T0PJwGDgdGBkPIz/qt34ghcfWNmg4AWzrnLi+13NTDb\nObe0hGNOBAaY2UJgNnAA/tGUz/cUTKrK+oxUyK1H04b1vC5FREQkrIRFOHHOTQgsXn0M/3TOAqCn\ncy490CQJaFV0HzOLB87Df8+TkjyOf6TkcWBfIB3/yExYrDlJ3ZFG1C5N6YiIiBQXFuEEwDk3AhhR\nyntXlrBtG1B/D8crDCaPB6vGYPorJ5XYAk3piIiIFFcTrtapljILUqn/jzXAIiIionDikZ2kkRit\nkRMREZHiFE48sisqlcZ1NXIiIiJSnMKJB3w+R0FsGkn1FU5ERESKUzjxwNr0TIjKpWVDTeuIiIgU\np3DigSVr/HeHbdtEIyciIiLFKZx44I9N/rvDtk/SyImIiEhxCiceWJXmHznp1EojJyIiIsUpnHhg\nfUYaFETRLqmh16WIiIiEHYUTD2zcnkpETlMiIszrUkRERMKOwokH0rNSqZOnKR0REZGSKJx4ICMv\njTi0GFZERKQkCice2OFSSYjUyImIiEhJFE48kBOZxj51FE5ERERKonDigbw6qTStr2kdERGRkiic\nVLEt27KhznZaJmjkREREpCQKJ1Vs6Vr/3WFb7aORExERkZIonFSxFev9d4fdP0kjJyIiIiVROKli\nq9L8IycdWiqciIiIlEThpIqt2eIfOTlw38YeVyIiIhKeFE6q2IbMVCx7H+Jio70uRUREJCwpnFSx\n9Kw0YnI1pSMiIlIahZMqtiU3lbo+XakjIiJSGoWTKra9II0GERo5ERERKY3CSRXLikilYYxGTkRE\nREqjcFLFcmNSaVpPIyciIiKlUTipQjm5+bjYv2gRr3AiIiJSGoWTKrRi3WYwR8tGmtYREREpjcJJ\nFVqx3n932P2aauRERESkNAonVej3Tf67wx7QXCMnIiIipVE4qUIpm/3hpHMbjZyIiIiURuGkCm3I\nTIPc+jROiPO6FBERkbClcFKFNu1IJWqXpnRERET2ROGkCm3JSSM2X1M6IiIie6JwUoUyC1Kpbxo5\nERER2ROFkyq0g1QSozVyIiIisicKJ1UoNyqNJnUVTkRERPZE4aSK+HyOgtg0mtXXtI6IiMieKJxU\nkZTUDIjKpVUjjZyIiIjsicJJFVm2zn/r+jaNNXIiIiKyJwonVeT3jYW3rtfIiYiIyJ4onFSR1en+\nkZNOrRRORERE9kThpIqs3ZoKBdG0aZbodSkiIiJhTeGkimzankpEdlMiIszrUkRERMKawkkV2Zyd\nRp18LYYVERHZG4WTKpKRl0o9tN5ERERkbxROqsgOl0ZCpMKJiIjI3kR5XUB1tWhVKms3Z5S5fVbU\nehrVOS6EFYmIiNQMCicVdPHIJ/m17ktl36EBtK3XJnQFiYiI1BAKJxU0pN9tLFnbp8ztoyMjufiU\n7iGsSEREpGZQOKmgU7u059Qu7b0uQ0REpMYJmwWxZnaTma0ys2wzm2VmR+yh7Vtm5jOzgsCvhV+/\nFWuXYGbDzWyDmeWY2TIzOyP0vQl/48aN87qEKqF+1izqZ82ifkppwiKcmFlf4AXgEeBwYCEwxcwa\nl7LLrUAS0Dzwa0tgCzChyDGjgWlAa+A/wIHAtcD60PSieqkt/1jUz5pF/axZ1E8pTbhM6wwARjrn\nxgKYWX/g38BVwLPFGzvntgPbC1+b2blAIjC6SLOrA9uOds4VBLatCUXxIiIiEjyej5wERji6Ad8U\nbnPOOfyjHseU8TBXAdOcc2uLbDsLmAmMMLNNZvabmd1nZp73WUREREoXDiMnjYFIILXY9lSgw952\nNrPmQC+gX7G39gNOBd4JvL8/8Ar+Pj9euZJFREQkVMIhnFTWFcBW4LNi2yPwB5zrAiMxv5hZS+Au\nSg8nsQBLly4NTaVhJDMzk/nz53tdRsipnzWL+lmzqJ81S5HvnbGVPZb5v297JzCtkwWc75z7vMj2\n0UCCc+68vey/AvjcOXdXse3TgVzn3OlFtp0BTALqOOfySzjWRcC7Fe+NiIhIrXexc+69yhzA85ET\n51yemc0DegCfA5iZBV4P29O+ZnYy0B54s4S3fwKSi23rAGwsKZgETAEuBlYDOWXrgYiIiOAfMWmL\n/3tppXg+cgJgZhfiv9KmP/Az/qt3LgA6OufSzWwQ0MI5d3mx/d4G2jvnji3hmC2BRcBY4CX8lxK/\nCbzonHs6hN0RERGRSvB85ATAOTchcE+Tx4BmwAKgp3MuPdAkCWhVdB8ziwfOw3/Pk5KOuc7MegJD\n8N83ZX3g9/+4NFlERETCR1iMnIiIiIgU0j0/REREJKwonIiIiEhYUTgJKM+DB6sjM3uk2EMSfWa2\nxOu6KsvMTjCzz81sfaBPZ5fQ5rHAwx+zzGyqme3vRa2Vtbe+FnkgZtGvL72qtyICd3H+2cy2mVmq\nmX1iZgeW0K5an9Oy9LOGnM/+ZrbQzDIDXzOKP3y1up/LQnvra004n8WZ2cBAPwYX217pc6pwQoUe\nPFhdLcK/4Dgp8HW8t+UERT38C6hvBP6xgMrM7gVuBq4DjgR24j+3MVVZZJDssa8Bk9n9HBe/nD7c\nnYD/6rqjgNOAaOBrM6tb2KCGnNO99jOgup/PtcC9QFf8jyn5FvjMzDpBjTmXhfbY14Dqfj7/FvgB\n/jr83y+Lbg/OOXXO1fovYBYwtMhrA9YB93hdWxD7+Agw3+s6QtxHH3B2sW0bgAFFXscD2cCFXtcb\ngr6+BXzsdW1B7mfjQF+Pr8nntJR+1rjzGejXX8CVNfVc7qGvNeZ8AvWB5fgfEfM/YHCR94JyTmv9\nyEmQHjxYXRwQmBL4w8zeMbNWe9+l+jKzdvh/Oil6brcBs6l557bQyYFpgmVmNsLMGnldUCUl4h8l\n2gI1+pzu1s8iasz5NLMIM+sHxAEzavC5/Edfi7xVU87ncGCic+7bohuDeU7D4j4nHqvUgwerkVn4\nn0O0HGgOPAp8b2YHO+d2elhXKCXh/w+/pHObVPXlhNxk4CNgFf47Jw8CvjSzYwKBu1oxMwNeBH50\nzhWuj6px57SUfkINOZ9mdjD+J8THAtuB85xzy83sGGreuSyxr4G3a8r57Ad0AbqX8HbQ/n0qnNQS\nzrmitxNeZGY/AynAhfiHG6Wac85NKPJysZn9BvwBnIx/6LW6GQF0Bo7zupAQK7GfNeh8LgMOAxLw\n3/l7rJmd6G1JIVNiX51zy2rC+TT/nddfBE5zzuWF8rNq/bQOsBkowL9IqahmwKaqL6dqOOcygRVA\ntVwZX0ab8K8fqlXntpBzbhX+v9/V7hyb2ctAb+Bk59zGIm/VqHO6h37+Q3U9n865fOfcn865X5xz\nD+BfQHkbNexcwh77WlLb6ng+uwFNgPlmlmdmecBJwG1mlot/hCQo57TWh5NA+it88CCw24MHZ5S2\nX3VnZvXx/6PY43+I1VngH/8mdj+38fivkKix57ZQ4Kecfahm5zjwDfsc4BTn3Jqi79Wkc7qnfpbS\nvlqezxJE4H8yfI05l3sQAdQp6Y1qej6nAYfgn9Y5LPA1F3gHOMw59ydBOqea1vEbDIw2/9ORCx88\nGIf/YYQ1gpk9B0zEP5WzL/BfIA/4v/buLVSqKo7j+PdXFpGZlZn6UJJZlD0UJlRGSAYFhkH10IVS\nErXoISoSslKju5GYlJegMhKLoAtW0g3swUoryKKypFBLDcvUPN66qP8e1hrbDufMuY3O9vj7wODM\nXvvyX7M85/xn7bVmvdLIuDpLUndSkqW8aYCkc4BNEbGG1AV5v6SfSKtNP0SaibWgAeF2Sq265scU\n0j3t9Xm/qaTesU6vEHqgSJpFml55JbBdUuUT2JaIqKwUftC3aWv1zG3dFdrzUdJYi1+AHqRV34cB\nl+VdDvq2rKhV167Snnl84j7fjyVpO7AxIr7Pm+rTpo2eklSWB+m7I1aTpjwtAYY0OqY61++V/B9k\nJ+mH52Xg1EbHVYd6DSNNwdxd9XihsM8DpOltO0i/CAY2Ou5615U0AO890i++v4CVwGygd6Pjbmcd\nm6vfbmBU1X4HdZu2Vs8u1J7P5dh35rp8AAzvSm3Zlrp2lfZsod6LKEwlrlebeuE/MzMzK5VDfsyJ\nmeCNPXgAAARcSURBVJmZlYuTEzMzMysVJydmZmZWKk5OzMzMrFScnJiZmVmpODkxMzOzUnFyYmZm\nZqXi5MTMzMxKxcmJmZmZlYqTEzOzVkhaJen2RsdhdqhwcmJmzZL0oqQ9eZG66rKZueyF/RzD6Hyd\n3ZJ2SdokaamkSXm10/1xvc31Pq+ZtY+TEzNrSZAWibxO0t5l3/Pz60krXB8IW4C+pNW0LwSeBUYB\nX0nqW+driVRvM2sgJydmVssyYA1wdWHb1aTEZFlxR0mXS1osabOkPyS9LWlAofwmSVslnVbYNkvS\ncklH1YghImJDRPwWESsiYi4wFDgGeKJwLkmaKGmlpB2Slkm6plA+LPfCjJD0taSdkpZIOrtSTlrh\nuWeht2ZyIY7ukp6X1CTpZ0nj2v42mll7ODkxs1qC9Ad7TGHbGGAuqZehqDswDRgMDAd2A2/uPVHE\nPGAh8LKkwyRdkc91Q0T81a6gIjYA84ErJVXiuBe4ERgPDAKmA/MkXVx1+BPAncAQYAPwlqTDgU+B\nO4AmoA/QD3iycNxdwBfAucAsYLak09sTt5m1TbdGB2BmpTcfeFzSyaQPNEOBa4FLijtFxBvF15LG\nAr9LGhQRy/PmW4GvgadJPTBTIuKrDsb1A9AD6CWpCZgIXBoRn+Xy1TkxuQVYXDjugYhYlGMcDawF\nroqI1yRtIffUNHO9hRExJz+fKulO0nvwYwfjN7MWODkxs5oi4g9J7wA3k3pLFkbEpv87LBJJA4EH\ngfOBE0mJTACnAMvzuf7MScv7wCcRMbUToVUCCGAgcDTwofYN7Ajgy2J1gKWFum2WtAI4qw3X+6bq\n9XrgpPYGbWatc3JiZm0xF3iG9Mf9thb2eQdYBYwFfiUlJ98BR1btNwzYBfST1D0itncwpkFAU0Rs\nLIxjGZGvXfR3B89f7d+q14FvjZvtF/7BMrO2eI+UZHQDPqgulHQCcAbwcER8FBErgF7N7DcUmACM\nBLYBMzsSjKSTSDOGKmNalpOSkP4RsbLqsa54KHBB4TzH57grt53+AQ7vSExmVj/uOTGzVkXEHkln\n5ufNTbXdDGwExktaD/QHHqMwLVdSD+AlYEZEvC9pHfC5pLcj4vUal5ekPqTE4jjSmJeJ+ZoTc0zb\nJD0JTM+DWz8GegIXAVvyYNyKyZI2Ab8Dj5AGxS7IZauBYyQNJ42N2RERO9v0JplZ3bjnxMzaJCK2\nRcS2FsqCNEj2PNLYjGnA3ZXi/O9TwFbgvnzMt/n5HEn9alz6WNKtmrWkGTXjSLeZBkfEb4UYJgEP\nAfeQekLeJd3mWVUMNZfPIM286Q2MjIhd+RxLgDnAq6TkZUJVHfapdo2YzawT1PyHIDOzriV/j8ki\n4PiIaGp0PGbWMvecmNmhpPq7WcyshJycmNmhxF3FZgcB39YxMzOzUnHPiZmZmZWKkxMzMzMrFScn\nZmZmVipOTszMzKxUnJyYmZlZqTg5MTMzs1JxcmJmZmal4uTEzMzMSuU/02lyoEQ9hsMAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff6eb4b3f28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(all_depth, all_l)\n",
    "plt.xlabel('Max Depth')\n",
    "plt.ylabel('Accuracy')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Trees = 128\n",
      "Trees trained.\n",
      "Total training time:  267.12233543395996\n"
     ]
    }
   ],
   "source": [
    "# train random forest for Kaggle submission\n",
    "n_trees=128\n",
    "n_sample=len(train_Xib)\n",
    "n_feature=-1 # -1 ==> sqrt(d)\n",
    "max_depth=80\n",
    "\n",
    "rf = random_forest(n_trees=n_trees, n_sample=n_sample, n_feature=n_feature, max_depth=max_depth)\n",
    "\n",
    "start = time.time()\n",
    "rf.train_parallel(train_Xib, train_y, verbose=False) # use all data\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# predict and save labels\n",
    "test_y = rf.predict(test_Xib)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame({'Category': test_y.astype(int)})\n",
    "df.index += 1\n",
    "df.index.rename('Id', inplace=True)\n",
    "df.to_csv('census_dist/census_rf_poly2_%s_binirize_nt%d_sample%d_m%d_depth%d.csv' \\\n",
    "          % (method, n_trees, n_sample, n_feature, max_depth))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Titanic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pclass</th>\n",
       "      <th>sex</th>\n",
       "      <th>age</th>\n",
       "      <th>sibsp</th>\n",
       "      <th>parch</th>\n",
       "      <th>fare</th>\n",
       "      <th>embarked</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3</td>\n",
       "      <td>male</td>\n",
       "      <td>?</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8.05</td>\n",
       "      <td>S</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>male</td>\n",
       "      <td>22</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>135.633</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>male</td>\n",
       "      <td>23</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>15.0458</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2</td>\n",
       "      <td>male</td>\n",
       "      <td>42</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>13</td>\n",
       "      <td>S</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3</td>\n",
       "      <td>male</td>\n",
       "      <td>20</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9.8458</td>\n",
       "      <td>S</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  pclass   sex age sibsp parch     fare embarked\n",
       "0      3  male   ?     0     0     8.05        S\n",
       "1      1  male  22     0     0  135.633        C\n",
       "2      2  male  23     0     0  15.0458        C\n",
       "3      2  male  42     0     0       13        S\n",
       "4      3  male  20     0     0   9.8458        S"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_X = pd.read_csv('./titanic_dist/titanic_training.csv')\n",
    "test_X = pd.read_csv('./titanic_dist/titanic_testing_data.csv')\n",
    "train_y = train_X.pop('survived').as_matrix()\n",
    "train_X.drop(['ticket', 'cabin'], inplace=True, axis=1)\n",
    "test_X.drop(['ticket', 'cabin'], inplace=True, axis=1)\n",
    "\n",
    "train_X.fillna(value='?', inplace=True)\n",
    "test_X.fillna(value='?', inplace=True)\n",
    "\n",
    "train_X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000, 7)\n",
      "(310, 7)\n"
     ]
    }
   ],
   "source": [
    "train_X = train_X.as_matrix()\n",
    "test_X = test_X.as_matrix()\n",
    "\n",
    "print(train_X.shape)\n",
    "print(test_X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000, 7)\n",
      "(310, 7)\n"
     ]
    }
   ],
   "source": [
    "# impute data\n",
    "train_Xi = imputer(train_X, cols=[1,6], method='max', missing='?')\n",
    "test_Xi = imputer(test_X, cols=[1,6], method='max', missing='?')\n",
    "method = 'median'\n",
    "train_Xi = imputer(train_Xi, cols=[0,2,3,4,5], method=method, missing='?')\n",
    "test_Xi = imputer(test_Xi, cols=[0,2,3,4,5], method=method, missing='?')\n",
    "all_Xi = np.vstack((train_Xi, test_Xi))\n",
    "\n",
    "print(train_Xi.shape)\n",
    "print(test_Xi.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# binarize categorical feature\n",
    "lb = LabelBinarizer()\n",
    "\n",
    "\n",
    "all_Xib = all_Xi[:, [0,2,3,4,5]]\n",
    "\n",
    "for j in [1, 6]:\n",
    "    b = lb.fit_transform(all_Xi[:, j])\n",
    "    all_Xib = np.hstack((all_Xib, b))\n",
    "    \n",
    "# polynomials for the numerical features\n",
    "# poly = PolynomialFeatures(2)\n",
    "# all_Xib = poly.fit_transform(all_Xib)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000, 9)\n",
      "(310, 9)\n"
     ]
    }
   ],
   "source": [
    "train_Xib = all_Xib[:train_Xi.shape[0], :]\n",
    "test_Xib = all_Xib[-test_Xi.shape[0]:, :]\n",
    "print(train_Xib.shape)\n",
    "print(test_Xib.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(800, 9)\n",
      "(200, 9)\n"
     ]
    }
   ],
   "source": [
    "# divide data into train/validation\n",
    "S = np.random.choice(range(len(train_y)), int(len(train_y)/5*4), replace=False)\n",
    "train_data = train_Xib[S, :]\n",
    "train_labels = train_y[S]\n",
    "S = np.setdiff1d(range(len(train_y)), S)\n",
    "validation_data = train_Xib[S, :]\n",
    "validation_labels = train_y[S]\n",
    "print(train_data.shape)\n",
    "print(validation_data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total training time:  0.14397740364074707\n",
      "Accuracy:  0.7925\n",
      "Accuracy:  0.79\n"
     ]
    }
   ],
   "source": [
    "dt = decision_tree(3)\n",
    "start = time.time()\n",
    "dt.train(train_data, train_labels)\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)\n",
    "\n",
    "l = dt.accuracy(train_data, train_labels)\n",
    "print('Accuracy: ', l)\n",
    "\n",
    "l = dt.accuracy(validation_data, validation_labels)\n",
    "print('Accuracy: ', l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dt.root.right.left.label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Trees = 128\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/working/lupolab/eason/anaconda3/envs/dsb2017/lib/python3.5/site-packages/ipykernel/__main__.py:31: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/working/lupolab/eason/anaconda3/envs/dsb2017/lib/python3.5/site-packages/ipykernel/__main__.py:31: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/working/lupolab/eason/anaconda3/envs/dsb2017/lib/python3.5/site-packages/ipykernel/__main__.py:31: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trees trained.\n",
      "Total training time:  3.306265354156494\n",
      "Accuracy:  0.92375\n",
      "Accuracy:  0.87\n"
     ]
    }
   ],
   "source": [
    "# train random forest\n",
    "n_trees=128\n",
    "n_sample=len(train_data)\n",
    "n_feature=-1 # -1 ==> sqrt(d)\n",
    "max_depth=30\n",
    "\n",
    "rf = random_forest(n_trees=n_trees, n_sample=n_sample, n_feature=n_feature, max_depth=max_depth)\n",
    "\n",
    "start = time.time()\n",
    "rf.train_parallel(train_data, train_labels, verbose=False) # use all data\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)\n",
    "\n",
    "l = dt.accuracy(train_data, train_labels)\n",
    "print('Accuracy: ', l)\n",
    "\n",
    "l = rf.accuracy(validation_data, validation_labels)\n",
    "print('Accuracy: ', l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Trees = 128\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/working/lupolab/eason/anaconda3/envs/dsb2017/lib/python3.5/site-packages/ipykernel/__main__.py:31: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trees trained.\n",
      "Total training time:  5.065195560455322\n"
     ]
    }
   ],
   "source": [
    "# train random forest\n",
    "n_trees=128\n",
    "n_sample=len(train_Xib)\n",
    "n_feature=-1 # -1 ==> sqrt(d)\n",
    "max_depth=80\n",
    "\n",
    "rf = random_forest(n_trees=n_trees, n_sample=n_sample, n_feature=n_feature, max_depth=max_depth)\n",
    "\n",
    "start = time.time()\n",
    "rf.train_parallel(train_Xib, train_y, verbose=False) # use all data\n",
    "end = time.time()\n",
    "print('Total training time: ', end-start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# predict and save labels\n",
    "test_y = rf.predict(test_Xib)\n",
    "df = pd.DataFrame({'Category': test_y.astype(int)})\n",
    "df.index += 1\n",
    "df.index.rename('Id', inplace=True)\n",
    "df.to_csv('titanic_dist/titanic_rf_max_binirize_nt%d_sample%d_m%d_depth%d.csv' % (n_trees, n_sample, n_feature, max_depth))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# writeup\n",
    "graph = pydot.Dot(graph_type='graph')\n",
    "graph.add_edge(pydot.Edge(\"sex, th=0.5 (male=1)\", \"pclass, th=2.5\"))\n",
    "graph.add_edge(pydot.Edge(\"sex, th=0.5 (male=1)\", \"age, th=9.5\"))\n",
    "graph.add_edge(pydot.Edge(\"pclass, th=2.5\", \"1 \"))\n",
    "graph.add_edge(pydot.Edge(\"pclass, th=2.5\", \"1  \"))\n",
    "graph.add_edge(pydot.Edge(\"age, th=9.5\", \"1\"))\n",
    "graph.add_edge(pydot.Edge(\"age, th=9.5\", \"0\"))\n",
    "# ok, we are set, let's save our graph into a file\n",
    "graph.write_png('graph.png')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
